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Abstract 

The goal of this paper is to quantitatively describe some statistical properties of higher- 
dimensional determinantal point processes with a primary focus on the nearest-neighbor distri- 
bution functions. Toward this end, we express these functions as determinants of x matrices 
and then extrapolate to — > cxd. This formulation allows for a quick and accurate numerical 
evaluation of these quantities for point processes in Euclidean spaces of dimension d. We also 
implement an algorithm due to Hough et. al. [1] for generating configurations of determinantal 
point processes in arbitrary Euclidean spaces, and we utilize this algorithm in conjunction with 
the aforementioned numerical results to characterize the statistical properties of what we call the 
Fermi-sphere point process for d = 1 to 4. This homogeneous, isotropic determinantal point process, 
discussed also in a companion paper [2] , is the high-dimensional generalization of the distribution 
of eigenvalues on the unit circle of a random matrix from the circular unitary ensemble (CUE). 
In addition to the nearest-neighbor probability distribution, we are able to calculate Voronoi cells 
and nearest-neighbor extrema statistics for the Fermi-sphere point process and discuss these as the 
dimension d is varied. The results in this paper accompany and complement analytical properties 
of higher-dimensional determinantal point processes developed in j2]. 
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I. INTRODUCTION 



Stochastic point processes (PPs) arise in several different areas of physics and math- 
ematics. For example, the classical statistical mechanics of an ensemble of interacting 
point particles is essentially the study of a random point process with the Gibbs mea- 
sure dfi{X) = P]\f{X)dX = exp[— j3V {X)]dX providing the joint probability measure for 
an A^-tuple of vectors X = (xi,...,Xiv) to be chosen. Moreover, some many-body prob- 
lems in quantum mechanics, as we will see, can be regarded as stochastic point processes, 
where quantum fluctuations are the source of randomness. With regard to mathematical 
applications, it has been well-documented [3] that the distribution of zeros of the Riemann 
zeta function on the critical line is well-represented by the distribution of eigenvalues of a 
random N x N Hermitian matrix from the Gaussian unitary ensemble (GUE) or circular 
unitary ensemble (CUE) in the limit oo. Nevertheless, it remains an open problem to 
devise efficient Monte Carlo routines aimed at sampling these processes in a computationally 
efficient way. 

In studies of the statistical mechanics of point-like particles one is usually interested in a 
handful of quantities such as n-particle correlation functions, the distributions of the spac- 
ings of particles, or the distributions of the sizes of cavities. Although these statistics involve 
only a small number of particles, it is not simple to extract them from knowledge of the joint 
probability density P/y. In general numerical techniques are required because analytical re- 
sults are rare. It is then of paramount importance to study point processes for which analytic 
results exist for at least some fundamental quantities. The quintessential example of such 
a process is the so-called Poisson PP, which is generated by placing points throughout the 
domain with a uniform probability distribution. Such a process is completely uncorrelated 
and homogeneous, meaning each of the ri-particle distribution functions is equal to p" , where 
p = N/V is the number density for the process. Configurations of points generated from 
this process are equivalent to classical systems of noninteracting particles or fully-penetrable 
spheres [Ij, and almost all statistical descriptors may be evaluated analytically. 

One nontrivial example of a family of processes which has been extensively studied is the 
class of determinantal PPs, introduced in 1975 by Macchi [5] with reference to fermionic 
statistics. Since their introduction, determinantal point processes have found applications 
in diverse contexts, including random matrix theory (RMT), number theory, and physics 
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(for a recent review, see [6]). However, most progress has been possible in the case of point 
processes on the hne and in the plane, where direct connections can be made with RMT [5] 
and completely integrable systems [7]. 

Similar connections have not yet been found, to the best of our knowledge, for higher di- 
mensional determinantal point processes, and numerical and analytical results in dimension 
d > 3 are missing altogether. In this paper and its companion [2J, we provide a generalization 
of these point processes to higher dimensions which we call Fermi-sphere point processes. 
While in [5] we have studied, mainly by way of exact analyses, statistical descriptors such 
as n-particle probability densities and nearest-neighbor functions for these point processes 
here we base most of our analysis on an efficient algorithm [1] for generating configurations 
from arbitrary determinantal point processes, and are therefore able to study other particle 
and void statistics related to nearest-neighbor distributions and Voronoi cells. 

In particular, after presenting in detail our implementation of an algorithm [1] to gener- 
ate configurations of homogenous, isotropic determinantal point processes, we study several 
statistical quantities thereof, including Voronoi cells statistics and distributions of minimum 
and maximum nearest neighbor distances (for which no analytical results exist). Addition- 
ally, the large-r behavior of the nearest-neighbor functions is computationally explored. We 
provide substantial evidence that the conditional probabilities Gp and Gy, defined below, 
are asymptotically linear, and we give estimates for their slopes as a function of dimension 
d between one and four. 

The plan of the paper is as follows. Section II provides a brief review of determinantal 
point processes and defines the statistical quantities used to characterize these systems. Of 
particular importance is the formulation of the probability distribution functions governing 
nearest-neighbor statistics as determinants of NxN matrices; the results are easily evaluated 
numerically. The terminology we develop is then applied to the statistical properties of 
known one- and two-dimensional determinantal point processes in Section III. Section IV 
discusses the implementation of an algorithm for generating determinantal point processes 
in any dimension d, and we combine the results from this algorithm and the numerics of 
Section II to characterize the so-called Fermi sphere point process for d = 1,2, 3, and 4. 
In Section V we provide an example of determinantal point-process on a curved space (a 
2-sphere) and our conclusions are collected in Section VI. 
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II. FORMALISM OF DETERMINANTAL POINT PROCESSES 



A. Definitions: n-particle correlation functions 

Consider point particles in a subset of d-dimensional Euclidean space E C M*^. It is 
convenient to introduce the Hilbert space structure given by square integrable functions on 
E; we will adopt Dirac's bra-ket notation for these functions. Unless otherwise specified, all 
integrals are intended to extend over E. A determinantal point process can be defined as a 
stochastic point process such that the joint probability distribution P/v of points is given 
as a determinant of a positive, bounded operator Ti. of rank A^: 

Pn{^i, xat) = det[f/'(xj, ^j)]i<i,j<N, (1) 

where (x, y) is the kernel of Ti. In this paper, we focus on the simple case in which the N 
nonzero eigenvalues of 7i are all 1; the more general case can be treated with minor changes 
[1]. We can write down the spectral decomposition of H as: 

N 

'>^ = T.\^nmi (2) 

n=l 

where {|(/>°)}^_-^ are the eigenvectors of the operator Ti.. The reason for the superscript 
on the basis vectors will be clarified momentarily. The correct normalization of the point 
process is obtained easily since [S]: 

/det|H(x..x,)].,„,„dx,...dx„ = iV!det(W). (3) 

where the last determinant is to be interpreted as the product of the non-zero eigenvalues 
of the operator Ti. Since these eigenvalues are all unity we obtain det(7i) = 1, which yields: 

PAr(xi, ...,XAr)ciXi ■ ■ ■ (ixAT = 1. (4) 

Notice that in terms of the basis we can also write: 

Piv(xi, ...,X7v) = — |det[0°(Xj)]i<ij<Ar| . (5) 

An easy proof is obtained by considering the square matrix <I>jj = 0°(xj) = ^Xj|0°). Then, 
|det[</.°(x,)]i<„<^|' = det($t)det(<l>) = det($t$) 



det 



N 



nl / l^j/ 



det[i7(x,,x,)], (6) 



which is the same as ([T]). 

Determinantal point processes are pecuhar in that one can actually write all the n- 
particle distribution functions explicitly. The n-particle probability density, denoted by 
p„(xi, . . . , x^), is proportional to the probability density of finding the first n particles 
in volume elements around the given positions (xi, . . . ,x„), irrespective of the remaining 
N — n particles. For a general determinantal point process this function takes the form: 

p„(xi, x„) = det[i^(xi, Xj)]i<ij<„ = n! P„,(xi, . . . , x„). (7) 

In particular, the single-particle probability density is: 

pi(xi) = i7(xi,xi). (8) 

This function is proportional to the probability density of finding a particle at xi, also known 
as the intensity of the point process. One can see that the normalization is: 

pi(x)rfx = tr(7^) = AT. (9) 

For translationally-invariant processes pi(x) = p, independent of x. We remark in passing 
that for a finite system translational invariance is defined in the sense of averaging the 
location of the origin over with periodic boundary conditions enforced. 
It is also possible to write the two-particle probability density explicitly: 

P2(xi,X2) = f/'(xi,Xi)if(x2,X2) - | i^'(xi , X2) | ^ , (10) 

which has the following normalization: 

j p2(Xi, X2)dXidX2 = N{N - 1). (11) 

In general the normalization for p„ is given by A^! / {N —n)\, or the number of ways of choosing 
an ordered subset of n points from a population of size A^. For a translationally-invariant 



and completely uncorrelated point process (10) simplifies according to p2 = p . 



We also introduce the n-particle correlation functions Qn, which are defined by: 

/ \ Pri(,Xi, . . . , X„j /1 o\ 

^n(Xl, . . . , X„) = — . (12) 

Since p„ = p" for a completely uncorrelated point process, it follows that deviations of 
from unity provide a measure of the correlations between points in a point process. Of 



particular interest is the pair correlation function, which for a translationally-invariant point 
process of intensity p can be written as: 

-f/'(xi,X2) 



5(2(Xi,X2) = = 1 



(13) 



P 

Closely related to the pair correlation function is the total correlation function, denoted 
by h\ it is derived from g2 via the equation: 



/i(x, y) = (72(x, y) - 1 = -p-' |//(x, y) 1^ , (14) 



where the second equality applies for all determinantal point processes by (13). Since 
5'2('") — ^ 1 as r — > cxD (r = ||x — y||) for translationally invariant systems without long-range 
order, it follows that h{r) ^0 in this limit, meaning that h is generally an function, and 
its Fourier transform is well-defined. 

Determinantal point processes are self-similar; integration of the n-particle probability 
distribution with respect to a point gives back the same functional form This property 
is desirable since it considerably simplifies the computation of many quantities. However, we 
note that even complete knowledge of all the n-particle probability distributions is not suffi- 
cient in practice to generate point processes from the given probability Pn- This notoriously 
difficult issue is known as the reconstruction problem in statistical mechanics [H |9l [TOl [11] . 
When in Section IV we discuss an explicit constructive algorithm to generate realizations 
of a given determinantal process, the reader should keep in mind that the ability to write 
down all the n-particle correlation functions Qn is not the reason why there exists such a 
constructive algorithm. 



B. Exact results for some statistical quantities 

We have seen that the determinantal form of the probability density function allows us to 
write down all n-particle correlation functions gn in a quick and simple manner. However, we 
can also express more interesting functions, such as the probability of having an empty region 
T) or the expected number of points in a given region, as properly constructed determinants 
of the operator Ti. This property has been used in random matrix theory to find the exact 
gap distribution of eigenvalues on the line in terms of solutions of a nonlinear differential 
equation [12]. The relevant formula is a special case of the result |6j that the generating 
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function of the distribution of the number points nj) in the region V is: 

= J2 ^(^^ = ^)'^" = + (-^ - l)xvnxv], (15) 

where xv is the characteristic function of V, I is the identity operator, and 2; G M. We 
will also denote P„ = P{nx) = n). Therefore, the probability that the region V is empty is 
obtained by taking the limit z — > in the previous formula. The result is: 

Po = det (I - xvHxv) • (16) 



Equation (16) may be written more explicitly. Consider the eigenvalues Aj of XvT^Xv- By 



the definition of the determinant, equation (16) takes the form: 



i=l 

where the product is over the non-zero eigenvalues of XvT^Xv only (of which there are A^, the 
number of particles). First notice that for the non-zero Aj we have Aj = Aj, where 
are the eigenvalues of Tixv'H. In fact one can show that the traces of all the powers of 
these two operators are the same using Xx> = Xvi 'H^ = 'H, and the cyclic property of the 
trace operation. This condition is sufficient for finite, and the limit ^ 00 can be taken 
afterward. The operator Tixv'H can now be written in a basis {(l)ri}n=i the matrix: 



Mi,{V) = / 0,(x)0,(x)rfx, (18) 



V 



and the determinant in (16) as: 

Po = det {5ij - M,j) . (19) 

We will be using this formula often in the following analysis. The probability Pq has a 
unique role in the study of various point processes [4], in particular when V = B{0; r), a ball 
of radius r (for translationally-invariant processes the position of the center of the ball is 
immaterial). In this context, Pq is called the void exclusion probability Ey{r) [H [T^ [T^ [T5]. 
and we will adopt this name and notation in this paper (in |2j we have studied this quantity 
in an appropriate scaling limit, when d ^ 00). 

However, there are statistical quantities of great importance which cannot be found with 
the above formalism. For example, one can examine the distribution of the maximum or 
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minimum nearest neighbor distances in a determinantal point process, or the "extremum 
statistics," and these quantities cannot be found easily by the above means. One could 
also explore the distribution of the Voronoi cell statistics or the percolation threshold for 
the PP. To determine these quantities we will have to rely on an explicit realization of a 
determinantal point process. The existence and the analysis of an algorithm to perform this 
task is a central topic of this paper. 

We introduce now some quantities which characterize a PP [U [131 [El US] . We start with 
the above expression Ey{r) for the probability of finding a spherical cavity of radius r in 
the point process. Analogously, one can define the probability of finding a spherical cavity 
of radius r centered on a point of the process, which we denote as Ep{r). Ep can be found 
in connection with Ey using the following construction. Consider the probability of finding 
no points in the spherical shell of inner radius e and outer radius r, which we call -Ey(r; e). 



This function can be obtained by either of the previous formulas (16) or (19). It is clear 
that Ey{r) = Ey{r; 0). It is also true that for sufficiently small e the probability of having 
two or more points in the sphere of radius e is negligibly small compared to the probability 
of having one particle. Hence, the probability Q{r; e) of finding no particles in the spherical 
shell -8(0; r) \ -8(0; e) conditioned on the presence of one point in a sphere of radius e and 
volume f (e) is: 

and by taking the limit e — >^ of this expression we find that: 

Ep{r) = \imn{r;e). (21) 

That Ep{0) = 1 can be seen from the following argument: set r = e + 0"^. Then, Ey{e + 
O'^; e) = 1 because the region is infinitesimal and hence empty with probability 1, and 
Ey{e; 0) ~ 1 — pv{e) since for sufficiently small e we have at most one point in the region. 
One line of algebra provides the result. 

Using this expression, we can derive an interesting and practical result for Ep. First, no- 



tice that Ey{r; e) contains the matrix Mij(r; e) defined by (18), which when e ^ becomes: 



M,,(r; e) ~ M,,(r) - t;(e)0,(O)0,(O). (22) 
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Moreover, if we assume that I — M is invertible, we can see that to first order in 5M: 



det(I - M + 5M) = exp[lndet(I -M + 6M)] = exp{tr[ln(I -M + 6M)]} 
~ exp{tr[ln(I - M)] + ti[5M{I - M)-^]} 
~ det(I - M){1 + tr[(5Af (I - M)-^]}. 



(23) 



From (23) we find the final result: 



Ep{r) = Ev{r) tr [A{1 - M)"^] 



(24) 



0, we have M 0, and £'p(0) = ti{A) 



where Aij = (f)i{0)(f)j{0)/ p. Notice that for r 
Ei l0i(O)lVp = H{0,0)/p = 1 as expected. 

These two primary functions can be used to define four other quantities of interest. Two 
are density functions: 

dEv{r) 



Hv{r) 
Hpir) 



dr 
dEp{r) 

dr 



(25) 
(26) 



which can be interpreted as the probability densities of finding the closest particle at distance 
r from a random point of the space or another random point of the process, respectively. 
The other two functions are conditional probabilities: 

Hv{r) 



Gv{r) 
Gp{r) 



ps{r)Ev{r) 
Hp{r) 



(27) 
(28) 



ps{r)Ep{r) ' 

which give the density of points around a spherical cavity centered respectively on a random 
point of the space or on a random point of the process. We note that s(r) is the surface 
area of the ci-dimensional sphere of radius r. We will study the behavior of these functions 
for some determinantal PPs in Sections III and IV of this paper. 



From the definitions in (25)- (28) in conjuction with (19) and (24), it is possible to express 



Hy, Hp,Gv, and Gp as numerically-solvable operations on N x N matrices. The results 
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are: 



Hv{r) 
Hp{r) 

Gvir) 
Gp{r) 



Evir) tr 



(I - M) 



dr 



Hv{r) tr [A{I-M)-'^] 
-Evir) tr 



A(I-M)-i^(I-M)-i 



ps[r 
Gvir) - 



tr 



(I - M)- 



d 



ps{r) 



dr 



dr 
^dM 
dr 



{Intr [Ail-M)-^]]. 



(29) 

(30) 
(31) 
(32) 



The form G'p(r) = Gvir) — G'(r) (which serves as a definition of G) in (32) is of particular 



interest. If the correction term G'(r) > for all r, positivity and monotonicity of Gp (which 
must be proven independently) are then sufficient to ensure that for appropriately large r, 
G'p(r) ~ Gvir) in scaling. Although we have been unable to develop analytic results for 



the large-r behavior of G, numerical results, which are provided later (see Fig. 10), suggest 
that G > and G ^ monotonically as r ^ oo for > 2, and G constant for d = 1. 
As both behaviors are subdominant with respect to the linear growth of Gy, we expect that 
Gp and Gv possess the same linear slope for sufficiently large r. 



An important point to address is the convergence of the results from (19) in the limit 



N —>■ oo. We expect that the calculations for finite but large provide an increasingly sharp 
approximation to the results from the N ^ oo limit. Figure [l] presents the calculation of Ep 
for a few values of A^ with (i = 1; it is clear that the numerical calculations quickly approach 
a fixed function for A^ > 40, and it is this function which we accept as the correct large- A^ 
limit. The results for higher-dimensional processes are similar, and we will assume that this 
convergence property holds throughout the remainder of the paper. 



C. Hyperuniformity of point processes 

Of particular significance in understanding the properties of determinantal point processes 
is the notion of hyperuniformity, also known as superhomogeneity. A hyperuniform point 
pattern is a system of points such that the variance o"^(i?) = (A^^) — (A^r)^ of the number 
of points Nn in a spherical window of radius R obeys: 

a\R) ~ i?^-^ (33) 
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FIG. 1: Convergence of the d = 1 numerical results using (19) for Ep with respect to increasing 
matrix size N. 

for large R [S]. This condition in turn implies that the structure factor S'(k) = 1 + ph{k.) 
has the following small-k behavior: 

lim S(k) = 0, (34) 
||k|Ho 

meaning that hyperuniform point patterns do not possess infinite-wavelength number fluc- 
tuations [8J. Examples of hyperuniform systems include all periodic point processes [8], 
certain aperiodic point processes [HI [IS], one-component plasmas [SI US], point processes as- 
sociated with a wide class of tilings of space [TTHIB], and certain disordered sphere packings 
[21 El UniiSO]. It has also been shown [2] that the Fermi sphere determinantal point process, 
described below, is hyperuniform. 



The condition in (34) suggests that for general translationally invariant nonperiodic sys- 
tems: 

S{k) r~.k" [k-^ 0) (35) 

for some a > 0. However, hyperuniform determinantal point processes may only exhibit 
certain scaling exponents a. One can see for a determinantal point process that: 

-p^{\H\'}{k\ 



S{k) 



(36) 



where ^ denotes the Fourier transform. Equation (|36|) therefore suggests that: 



{l-k") (k-^O). 



(37) 
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Taking the inverse Fourier transform of (37) gives the following large-r scaling of \H{r 



\H(r) 



pn 



d/2 



2" r(^^±^^) 



OO 



(38) 



The negative coefficient and the negative argument of the gamma function in (38) are cru- 



cially important. Since \H{r)\ > for all r, it must be true r(— a/2) < 0, and this condition 
restricts the possible values of the scaling exponent a. Namely, the behavior of the gamma 
function requires that a fall into one of the intervals (0,2], [4,6], [8,10], and so forth. We 
remark that the integer-valued endpoints of these intervals are indeed valid choices for a 
and imply that |-ff(r)| ~ for sufficiently large r. These values of a are therefore types 
of "limiting values" that overcome the otherwise dominant r~^°''^'^^ asymptotic scaling of 
\H{r)f. We provide an example of a determinantal point process with the critical scaling 
a = 2 in Section III.C; the resulting large-r behavior for H{r) is seen to be exponential. 

III. PROPERTIES OF KNOWN DETERMINANTAL POINT PROCESSES 
A. One-dimensional processes 

By far, the most widely studied examples of determinantal point processes are in one- 
dimension. In fact, the connection to (RMT) led others to explore the statistical properties 
of these systems even prior to the formal introduction of determinantal point processes. To 
make this connection explicit, consider an x random Gaussian Hermitian matrix, i.e., a 
matrix whose elements are independent random numbers distributed according to a normal 
distribution. This class of matrices defines the Gaussian unitary ensemble or GUE. It is 
possible to see [3] that the distribution induced on the eigenvalues of these random matrices 
is: 



PAr(Ai,...,A7v) = ^ TT |Ai - Ajpexp 



(39) 



where Zm is an appropriate normalization constant. By a standard identity for the Vander- 
monde determinant: 



niA. 

Kj 



A, 



det(A^ 



l<i,n<N\ , 



and by combining the rows of the matrix A" appropriately we find: 



^<i,j<N\ ; 



(40) 



(41) 



i<j 
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where the functions Hn{x) are the Hermite orthogonal polynomials normalized such that 

2 

the coefficient of the highest power of Hn is unity. Taking into account the weight , 
we can write in agreement with (j5]): 

Pn = ^1 det[0j(Ai)]i<i<j<Ar|^ (42) 

where the orthonormal basis set is: 

= -^H„{x) exp[-xV2]; (43) 



Zn is a normalization factor. Therefore, this distribution is equivalent to the the one induced 
by a system of non-interacting, spinless fermions in a harmonic potential. We note with- 
out proof that the other canonical random matrix ensembles (GOE and GSE) can also be 
expressed as determinantal point processes by introducing an internal vector index for the 
basis functions [3, 6J. 

Another prominent example of a ci = 1 determinantal point process is given by the unitary 
matrices distributed according to the invariant Haar measure; the resulting class is termed 
the circular unitary ensemble or CUE [21j. The eigenvalues of these matrices can be written 
in the form Xj = e*^^ with 9j G [0, 27i] Vj G N; they are distributed according to ^ with the 
basis: 

MO) = ^exp[tne]. (44) 
V 2tt 



Notice that the eigenvalues represent the positions of free fermions on a circle, where the 
Fermi sphere has been ffiled continuously from momentum to — 1. 

Another possible one-dimensional process is obtained by changing the exponent x"^ in 



(39) to an arbitrary polynomial. This generalization has interesting connections to the 
combinatorics of Feynman diagrams and to random polygonizations of surfaces [22]. For 
other examples of one- dimensional determinantal point processes we refer the reader to [6] . 



B. Exact results in one dimension 



For historical reasons, the most-studied descriptor of determinantal point processes is 
the gap distribution function, which represents the probability density of finding a chord of 
length s separating two points in the system for d = 1; we denote this function by p{s). 
For canonical ensembles of random matrices exact solutions for p{s) have been written in 
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terms of solutions of well-known nonlinear differential equations |3]. We start with the 
following observation: after an appropriate rescaling of the eigenvalues, the gap distribution 
of eigenvalues of a random matrix is a universal function, depending only on the "nature" 
of the ensemble (unitary, orthogonal or symplectic) which defines the small-r behavior of 
g2- For example, the two ensembles GUE and CUE defined above will have the same 
gap distribution in the limit N ^ oo. In the case of the GUE the limit is taken for the 
eigenvalues: 

A^ = ^ + -^^yi^ (45) 

where z is in the "bulk" of the distribution {\z\ < \/2N — e for large). One can prove 
that all the eigenvalues of a large random matrix will fall in an interval of size 2\/2N with 
probability 1 in the large- limit. After this rescaling, the kernel H converges to the "sine 
kernel" in the large- A^ limit [El EH]: 

Hn{Xi,X2) H{yi,y2) = — . (46) 

N~.oc ■n-{yi-y2) 

From this result one can find the n-particle correlation functions. In particular one finds for 

92- 

Applying this procedure to the CUE leads to the very same kernel; for a wider class of 
examples relevant to physics see [23]. Convergence of the kernel implies weak convergence 
of all the n-particle correlation functions to universal distributions. These distributions are 
defined by the sine kernel, one of a small family of kernels which appear to be universal 
[I2I [23] in controlling large- A^ limits of various statistical quantities of apparently different 
distributions. The study of the analytic properties of the kernels in this family yields a 
complete solution for the Janossy probabilities and edge distributions in one-dimensional 
systems. 

Once the limiting kernel is identified, a solution for the gap distribution p{s) still requires 
a detailed mathematical analysis [T^]. An approximate form for p{s), known as Wigner's 
surmise, was suggested by Wigner in 1951: 



2 



32s 

p[s) = — ^ exp 



TT 



TT 



(48) 



and it is an extremely good fit for numerical data. However, our primary focus in this work 
is on the asymptotic behavior of the conditional probability Gy, and we therefore look for 
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an exact solution for this function. First, we note without proof [25] that Ev{s) for c? = 1 
may be expressed in terms of a Painleve V transcendent. Namely, let a{s) be a solution of 
the nonlinear equation: 



{say + 4(sa' - a) [sa' - a + (a')'] = 0, 
subject to the boundary condition: 



(49) 



a[s) ~ — 

as s ^ 0. We may then write Ev{s) in the form: 

Ev{s) = exp 



(50) 









[Jo 


t 





(51) 



We recall that Ev{s) may also be expressed in terms of Gy via the following relation: 



Ey{s) = exp 



—2 / GY{x)dx 





(52) 



By making a change of variables and comparing (|51|) and (|52|), we conclude that: 

Gvis) : 



5-(27rs) 



2s 



(53) 



Equation (53) allows us to develop small- and large-s expansions of Gy in terms of the 
equivalent expansions for a. 

To describe the small-s behavior of Gy, we substitute an expansion of the form: 



S / S ^ 2 
TT VtT 



N 



n=3 



(54) 



into (|49|) and solve order-by-order for the coefficients 6„. Upon converting the solution to a 

167r2 647r^^ 



result for Gy using (53), we obtain: 

Stt^ 



Gy(s) = l + 2s + 4s" + 



9 



16 



207r^ 
~9~ 



32 



+ 



225 



1127r2 4487r^\ « ,^,7, 
+ ( 64 — + \ + Ois'). 



(55) 



9 675 

The derivation of the large-s expansion is similar. We choose an expansion of the form: 



N 



a{s) = h^s^ + hiS + 62 + ^ Ks 



2-n 



(56) 



n=3 



16 




Exact 

Small-r expansion 
Asymptotic expansion 



I I I I I I I I I I I I I I I I 
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 



FIG. 2: Comparison of the exact form of Gy for the d = \ determinantal point process with the 



small- and large-r expansions in (55) and (57) 



and substitute this equation into (49). After converting the result to an asymptotic series 



for Gy with (53), we obtain: 

. TXH 1 



G 
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6575 



1080091 16483607 



2 8s 327r2s3 QAti^s^ 2hQ'K^s'^ 10247r8s9 81927riOsii AmQ-n^^s^^ 



-15\ 



(57) 

By looking at Figure [2| one can see that the expansions are quite good for the ranges in s 



where they are valid. Equations (49), (55), and (57) constitute the solution to our problem 



Although it is natural to ask if there is a corresponding nonlinear differential equation that 
characterizes Gy in higher dimensions, we are not aware of any work in this direction, and 
this issue remains an open problem. 



C. Two-dimensional processes 



There are a few examples of determinantal point processes in two dimensions. The 
seminal example is provided by the complex eigenvalues of random non-Hermitian matrices 
[26t [27] . The kernel of such a determinantal point process is given by: 

N-l 



Hn{z,w) 



exp 



1 



\w\ 



E 

fc=0 



{zwy 



(58) 
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where is the rank of the matrix and z,w & C Incidentally, (58) can be related to the 
distribution of polarized electrons in a perpendicular magnetic field, filling the N lowest 



Landau levels. In the limit — oo (|58|) becomes 

1 



H(z,w) 



TT 



exp 



--(1^1^ + \w\^ -2zw) 



(59) 



which is a homogeneous and isotropic process (p = H{z,z) = ^) in C. It is instructive to 
examine the pair correlation function, which after some algebra can be written as: 



92(^1, Z2) = 1- exp {-\zi - • 



(60) 



From this expression one finds that the correlation between two points decays like a Gaussian 



with respect the distance separating the points. Letting r 
associated structure factor of the system as: 



S{k) = 1 — exp 



' 4 



1 2:1 — 2:2!, we may write the 



(61) 



O (A;^) (k ^ 0). 



(62) 



which has the following small-A; behavior: 

S{k) - ^ 

We see that the determinantal point process generated by the Ginibre ensemble is hyper- 
uniform with an exponential scaling a = 2 for small k, corresponding to an endpoint of one 
of the "allowed" intervals for determinantal PPs; the large-r behavior of the kernel H{r) is 
exponential {H{r) = exp[— r]). 

Other ensembles of two-dimensional determinantal point processes can be found in 
simple systems. For example, the n zeros of an analytic Gaussian random function 
akz'' also form a determinantal point process on the open unit disk [28l |29] . 
The limiting kernel governing these zeros is called the Bargmann kernel: 



H{Zi,Z2) 



[1 - ZIZ2] 



(63) 



and is inherently different from (59). 



IV. AN ALGORITHM FOR GENERATING DETERMINANTAL POINT PRO- 
CESSES 

We are able to write down an algorithm, which we call the HKPV algorithm, after 
[1], to generate determinantal point processes due to the geometric interpretation of the 
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determinant in as the volume of the simplex built with the vectors Vj = {|0j)}i<j<Af- 
In the original paper [T] this algorithm is sketched and then proved to produce the correct 
distribution function pi^j. The algorithm is extremely powerful and versatile and we believe 
it is important to provide as many details as possible about it and its implement ant ion 
(which has not been done before, to our knowledge). Therefore we dedicate the present 
section to provide a complete description of the HKPV algorithm and enough details (with 
some tricks) for its efficient implementation. 

Set ifiv = H, the kernel of the determinantal point process. Pick a point distributed 
with probability: 

p^(x) = i7;v(x,x)/iV. (64) 
With this point build the new operator A^^i, defined by: 

An.i = Hm\^n){^n\Hn. (65) 

This operator has with probability 1 a single nonzero eigenvalue and — 1 null eigenvalues. 
When expressed as a matrix in the basis {|0°)}i<n<Af, ^Af-i takes the form: 



(A^_i),^^. = mN^ii^N). (66) 



Consider the — 1 null eigenvectors of Atv-i; we will denote them as and call 

\4>n) the only eigenvector with a nonzero eigenvalue. The null eigenvectors can be found 
easily by means of a fast routine based on singular value decomposition (SVD), but we will 
see one that can proceed without it. 
Next, build the new operator Hn^i. 

Hn-i = I'/'nX^nl ) H^. (67) 



,n=l 

N 



To simplify the computation, notice that by completeness of the basis in the 

eigenspace of Hn: 

Hn [Y. \^n){4>l\ + \4>\){4>\\\ Hn = H^, (68) 



,ra=l 



and since |0]v) only eigenvector orthogonal to the null space: 

= tr(A^_i)|0jv>(0]vl- (69) 
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From this equation we conclude: 




tr(Ajv-i) 



1 



(70) 



Once i^AT-i is obtained, we repeat the procedure with H^^i, generating the point 

^Ar_i from the probabihty distribution: 



and the operators An_2, i^Ar_2- As the number of iterations increases, we constantly reduce 
the rank of the operators by one: tr(if7v) = N, tr(/77v-i) = — 1, etc. Therefore, after we 
have placed the last point ^i, we are left with the an operator of rank 0, and the algorithm 
stops. Reference [1] shows that the A^-tuples (^i,. . . ,^n) are distributed according to the 



The whole procedure requires O {N"^) steps for every realization, which is equal to the 
number of function evaluations necessary to create the matrices A. Therefore, the algorithm 
is computationally quite light. The only subroutine that requires some work is the extraction 
of the random points from the probability distributions Pn(x). For d = 1 one can use a 
numerically-implemented inverse CDF technique [30], and the computational cost of this 
procedure is independent of A^. For > 2 if the distributions are not very peaked, a 
rejection algorithm is sufficient. The rejection algorithm works by sampling points from a 
uniform distribution on the domain. A tolerance value near the maximum of the probability 
density of the point process is set, and the point is accepted if a uniform random number 
chosen between and the tolerance value is less than the probability density at that point. 
Otherwise, the point is rejected, and the process repeats. Unfortunately, it is difficult to 
estimate the computational cost of this algorithm as a function of the number of particles 



A. Numerical results in one dimension 

We have implemented the algorithm described above to study a determinantal point 
process on the circle x G [0,27r], where (pni^) = exp [ma;] / v^27r are the orthonormal 
functions with n = 0, ±1, ±2, . . . , ±N/2. This ensemble, as mentioned above, is equivalent 
to the one generated by the eigenvalues of unitary random matrices chosen according to the 



PAr_i(x) =i7^_l(x,x)/(Ar-l) 



(71) 



distribution ([T]). 



N. 
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Haar measure. Eigenvalues of matrices from the CUE can be generated easily by means of 
a fast SVD algorithm [21]; however, plenty of exact results exist. Therefore, we study this 
one-dimensional determinantal process as a test both for the performance of our algorithm 
and for the convergence of the results to the N ^ oo limit. 

We have implemented the algorithm in both Python and C++, noticing little difference 
in the speed of execution, and run it on a regular desktop computer. As mentioned above, 
the algorithm runs polynomially in with the sampling of the distribution pn limiting the 
computational speed. One point, however, which requires attention is the loss of precision 
of the computation. Due to the fixed precision of the computer calculations, the matrix Hn 
ceases to have exactly integral trace, diminishing the reliability of the results. Typically, one 
observes deviations in the 5th decimal place after ^ 50 particles have been placed. We have 
devised an 'error correction' procedure in which the numerical matrix Hn is projected onto 
the closest Hermitian matrix Hn which has eigenvalues 1 or only. This projection corrects 
for a great part, but not all, of the error; however, the algorithm is slowed by this modifi- 
cation. The number of particles in each configuration can therefore be pushed to ~ 100, 
regardless of the dimension. We have been able to generate between 75000 and 100000 con- 
figurations of points in each dimension. In general, the error-correction procedure generates 
more reliable statistics for a given value of A^ compared to the uncorrected algorithm, and we 
therefore expect that any residual error not captured by the matrix projection is minimal. 

As a preliminary check for the HKPV Algorithm, we have calculated the pair correlation 



function g2 and compared the results to the exact expression in ( 75 ) below. The comparison 
is quite favorable and suggests that the point configurations are being generated correctly 
by the algorithm. We mention a few characteristics of g2 which arise from the determinantal 
nature of the point process. First, the system is strongly correlated for a significant range 
in r, and g2{r) — ^ as r — » 0. This correlation hole [HI |32l [331 El] is indicative of a strong 
effective repulsion in the system, especially for small point separations. In other words, the 
points tend to remain relatively separated from each other as they are distributed through 
space. Second, (72 < 1 for all r, meaning that it is always negatively correlated; again, this 
quality is indicative of repulsive point processes, which are characterized by a reduction of 
the probability density near each of the coordination shells in the system. We show in a 
separate paper [2] that at fixed number density, g2 approaches an effective pair correlation 
function (^) ~ ^l^" — D) as d ^ 00, suggesting that the points achieve an increasingly 
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FIG. 3: Comparison of the exact expression (75) for g2{r) with the results from the HKPV Algo- 
rithm for d = 1, p = 1- The results from the simulation are obtained using 75000 configurations of 
45 particles. 

strong effective hard core D (~ \/d) as the dimension of the system increases. At fixed mean 
nearest-neighbor separation this observation implies that g2{r) — > ^2(^) ~ ^ ^"^^ all r > as 
d oo (as g2{0) = for any d), implying that the points become completely uncorrelated 
in this limit. We will show momentarily that the latter limit is difficult to interpret due to 
the dimensional dependence of the density p. 

Figure |4] presents the results for the gap distribution function p(r) using both the HKPV 
Algorithm and a numerical calculation based on the determinant in ( [l9| ). As with the calcu- 
lation of g2, the comparison between the numerical results and the simulation is favorable. 
This curve, as expected, has the same form as the one reported in the random matrix liter- 
ature [3j and scales with as r — > 0. We stress, however, that this function represents the 
distribution of gaps between points on the line and does not discriminate between gaps to 
the left and to the right of a point. The random matrix literature oftentimes describes this 
quantity as a "nearest-neighbor" distribution, which it is not. As mentioned in the discus- 



sion following (25) and (26), the void and particle nearest-neighbor distribution functions 



are given by the functions Hy and Hp, respectively, and require that distance measurements 
be made both to the left and to the right of a point; the numerical and simulation results 
for these functions are also given in Figure |4} 
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FIG. 4: Comparison of numerical and simulation results with d = 1, p = 1 for left: the gap 
distribution function p{r), center: Hp{r) and right: Hv{r). 

The function Hy is clearly different from p. Hp has a similar shape to the gap distribution 
function; however, Hp peaks more sharply around r ^ 0.725 while p has a less intense 
peak near r ^ 1. This observation is justified from a numerical standpoint since point 
separation measurements are made in both directions from a given reference point with only 
the minimum separation contributing to the final histogram of Hp. In contrast, every gap 
in the point process is used for constructing the histogram of p. As a result, we expect the 
first moment of Hp to be less than that of p, and this result is exactly what we observe in 
Figure |4} 

The form of Hy may at first seem confusing in the context of our discussion above 
concerning the inherent repulsion of the determinantal point process. Unlike Hp and p, the 
void nearest-neighbor function Hy has a nonzero value at the origin and is monotonically 
decreasing with respect to r. To understand this behavior, it is useful to examine the 
behavior of the corresponding Gy and Gp functions, which are plotted in Figure [5] We 
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FIG. 5: Numerical results using (19) for Gp and Gy with d = 1, p = 1. Also included are 



representative simulation results and estimated errors from the HKPV Algorithm under the same 
conditions as in Figure |3j 



recall from (27) and (28) that Gy and Gp are related to conditional probabilities which 



describe, given a region of radius r empty of points (other than at the center for Gp), the 
probability of finding the nearest-neighbor point in a spherical shell of volume s{r)dr, where 
s(r) is the surface area of a d-dimensional sphere of radius r. Of particular relevance to the 
behavior of Hy is the fact that Gv{0) = 1 and s(0) = 2 for = 1. Therefore, the dominant 
factor controlling the small-r behavior of Hy is the spherical surface area s(r) p]. Since 



s(0) is nonzero for d = 1, it follows from (27) that Hy{0) is nonzero in contrast to Hp{0). 



The behavior of both Gp and Gy is of particular interest in this paper. We conjecture 
that both functions are linear for sufficiently large r in any dimension. We show elsewhere 
[2] that as r ^ 0, Gp{r) ~ ^{dy + O (r^) and Gv{r) ~ 1 + C (r^), where ^{d) is a 
dimensionally-dependent constant (for d = 1 this is evident in Figure [sj. Additionally, we 
believe that Gy and Gp obtain the same slope in the large-r limit, and we will provide 



further commentary on this notion momentarily (see Fig, 10). It is clear from Figure |5] that 
the results from the simulations are in agreement with the numerical results for a wide but 
limited range of r and they begin to deviate respectably for r sufficiently large. This is due 
to the fast decay of both Hp/y and Ep/y to zero, therefore giving very small statistics (and 
a large degree of uncertainty) at these values of r. This said, the numerical results are clear 
and provide strong support for our claims above. 
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B. Fermi sphere deter minantal point process for d > 2 

1. Definition of the Fermi sphere point process 

Here we study the determinantal point process of free fermions on a torus, filling a Fermi 
sphere. A detailed description of this process in any dimension may be found in an accom- 
panying paper [2] . We consider this example because it is the straightforward generalization 
of the one-dimensional CUE process described above. However, sampling of this ensemble 
cannot be accomplished with methods other than the algorithm introduced above; this lim- 
itation is in contrast to the two examples from Section III.C, where the ensemble may be 
generated from zeros of appropriate random complex functions. Nevertheless, it is difficult 
to construct another procedure that can be generalized to higher dimensions since zeros of 
complex functions and random matrices are naturally constrained to d <2. 

We consider the determinantal point process obtained by "filling the Fermi sphere" in 
a d-dimensional torus, i.e., x G [0,27?]^^; our choice of the box size is for convenience and 
without loss of generality. We therefore consider all functions of the form: 

0n=(^j exp[?(n,x)] (72) 

with 

||nf < 4(iV), (73) 

where k']p{N) is implicitly defined by the total number of states contained in the reciprocal- 
space sphere. This process is translationally-invariant for any A^, both finite or infinite, and 
isotropic in the limit ^ cxd; it possesses the symmetry group of the boundary of the set 



(73), a dihedral group which approximates SO{d) very well for A^ sufficiently large. The 
pair correlation function can be easily calculated for any N < oo, and it is well-defined in 
the thermodynamic limit: 

92{x) = 1 - ^ 5^ 5^exp[z(n - n',x)], (74) 

n n' 



where n and x are c?-dimensional real vectors, and the sums extend over the set (73), which 
contains A^ points. In the limit N oo the sums become integrals over a sphere of radius 
kp = 2^/^T[pT{l + where p = N/{2tiY is the number density. The resulting pair 
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correlation function is given by: 

92{r) = 1 



2'^[r(i + d/2)]^ 



[Jd/2{kFr)f 



(75) 



where Jd/2 is the Bessel function of order d/2 (cf. [2]). This pair correlation function is 
clearly different from (60) for d = 2; the two are therefore not equivalent, even in the 
thermodynamic limit. One can also find the limiting kernel: 

-(2'i/2)r(l + c//2) 



{kF\\^-y\\r/^ 



Jd/2ikF\\^-y\ 



(76) 



which is also different from (59) and (63) for d = 2. 



Figure |6] shows configurations of points generated for the d = 2 and d = 3 Fermi sphere 
point process alongside corresponding configurations for the Poisson point process in these 
dimensions. The repulsive nature of the determinantal point process is immediately apparent 
from these figures; note especially that the Fermi sphere point process discourages clustering 
of the points in space. In contrast, clustering is not prohibited in the Poisson point process, 
and small two- and three-particle clusters are easily identified. Of particular interest is that 
the Fermi sphere point process distributes the points more evenly through space due to 
the effective repulsion in the system. This characteristic reflects the hyperuniformity of the 
point process [8J, and we will have more to say about this property momentarily. 



2. Calculation of §2 and nearest-neighbor functions for d >2 

Figure [7] shows the numerical and simulation results for the pair correlation function g2 
with d = 2] a comparison of the results provides strong evidence that the HKPV Algorithm 
correctly generates configurations of points for the Fermi sphere point process even in higher 
dimensions. Note that the d = 2 correlations are significantly diminshed with respect to the 
form of g2 for d = 1; this behavior is in accordance with a type of decorrelation principle 
[351 136] for the system. Namely, we expect that as the dimension of the system increases, 
unconstrained correlations in the system diminish. We also remark that all higher-order 
correlation functions Qn can be written in terms of the pair correlation function g2 for any 
determinantal point process. We prove this claim in an accompanying paper |i2j. It is 
therefore clear that the HKPV Algorithm is a powerful method by which one can study 
determinantal point processes in higher dimensions. 
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FIG. 6: Upper: A d = 2 configuration of = 109 points distributed according to left: a Fermi 
spliere determinantal point process and right: a Poisson point process. Lower: A d = 3 configura- 
tion of = 81 points distributed according to left: the Fermi sphere determinantal point process 
and right: a Poisson point process. All configurations have /? = 1. 

Figure |8] contain results for the nearest-neighbor particle and void density functions Hp 
and Hy for d = 2,3, and 4. In all cases the numerical results coincide with the simulation 
results. We do note that for d = 3 and = 4 we have implemented the error-correction 
procedure described in Section IV.A to increase the reliability of the simulation results as 
well as the particle numbers. As mentioned above, running the algorithm without error 
correction generally results in a loss of precision in the trace of the kernel matrix H during 
computation; the error introduced by this loss of precision as measured by deviation from 
the "exact" numerical results increases with respect to increasing particle number, and we 
notice that the errors are more acute for d = 3 and d = 4. Although some error still remains 
in the results even after projecting the matrix H onto the nearest Hermitian projection 
matrix, the results in these figures leave us confident that the computations are reliable. 

In contrast to the d = 1 process, Hy for d = 2,3,4 approaches as r ^ 0; for d = 3 
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FIG. 7: Comparison of the exact expression (75) for g2{r) with the results from the HKPV Algo- 



rithm for d = 2, p = 1- The results from the simulation are obtained using 75000 configurations of 
45 particles. 
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FIG. 8: Comparison of numerical and simulation results for left: Hp with d = 2,3, and 4 at number 
density p = 1 and right: for Hy with d = 2,3, and 4 at number density p = 1. 



and 4, Hy and Hp in fact possess very similar overall shapes. The small-r behavior of ifv 



in these cases is due to the behavior of s(r) for d > 2; namely, s{r) 



d-1 



for all d, and 



for > 2 we observe that s(0) = as opposed to the d = 1 case, where s(0) = 2. We have 
already shown with generality that Gv{r) — > 1 as r — > 0, a result which may be observed 
in Figures (qI One can see from these figures that Gy^r) — >■ 1 as r — in each dimension, 



28 



reinforcing the dominance of s(r) in the small-r behavior of Hv{r). 

For d = 2, the shape of Hy resembles the corresponding curve for a Poisson point process; 
nevertheless, these two processes are inherently different. We may easily see the deviation 
between the two processes by noticing that Hp and Hy do not coincide for any dimension 
and that Gp and Gy both increase linearly for sufficiently large r. The latter observation 
implies that Hv{r) ^ s{r)Ey{r) for large r, which is the case for the Poisson point process. 
However, we show elsewhere [2] that Ey for the Fermi sphere point process in dimension d 
(finite) behaves similar to the corresponding function for a Poisson point process except in 
dimension d + 1. Further justification for this claim is also developed later is this paper. 

With regard to iJp, we remark that in each dimension Hp{0) = in agreement with the 
repulsive nature of the point process. However, it is worthwhile to note that, in light of 
the connection to noninteracting fermions described above, we can associate this repulsion 
with a type of Pauli exclusion principle, which for noninteracting fermions is purely quantum 
mechanical in nature and arises solely from the constraint of antisymmetry of the A^-particle 
wavefunction. The determinantal form of the wavefunction is the manifestation of this 
antisymmetry in any dimension, thereby providing some physical insight into the strong 
small-r correlations for this determinantal point process. We stress that in the case of 
noninteracting fermions the repulsion does not arise from any true interaction among the 
particles and is purely a consequence of the aforementioned antisymmetry. 

We show in an accompanying paper [2] that for any d, Hp ~ r'^^^ for small r, and we 
observe this behavior in our results. It is also true [2] that Hy ~ r'^-^^ Ey ~ 1 — x{d)r'^, and 
Ep ~ 1 — x{d)r'^~^'^ as r ^ 0, where x{d) and >c{d) are dimensionally- dependent constants. 
These properties imply that Gp ^ r"^ and Gy ~ 1 for small r as with the d = 1 case. Figure [o] 
shows these trends in greater detail. With regard to the large-r behavior of Gp and Gy, the 
linearity of both curves apparently holds in each dimension. A surprising detail, however, is 
that Gp and Gy appear to converge with respect to increasing dimension. To understand 



this observation, we recall from (32) that Gp = Gy — G; Figure 10 provides plots of G{r) for 
d = 1, 2, 3, and 4. It is clear from these curves that in each dimension G for large r is positive 
and scales more slowly that r in each dimension. We therefore expect that the large-r slope 



of Gp is equal to the asymptotic slope of Gy according to (32). Since numerical results for 
Gy are more easily and more accurately obtained, we assume this asymptotic convergence 
and provide results for the asymptotic slope of Gy below. Table [T] collects our calculations 
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FIG. 9: Numerical results using (19) for Gp and Gv with (from left to right) d = 2, 3, 4 and for all 
p = 1. Also included are representative simulation results and estimated errors from the HKPV 
Algorithm. 




FIG. 10: Plots of G{r) = Gv{r) - Gp{r) for d = 1, 2, 3 and 4. 
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for the slope of Gy in each dimension for large r. The slopes are calculated by fitting the 
large-r portion of each quantity to a function of the form: 

F {x) = aox + ai + ^ . (77) 

It has been conjectured in j2] that as the dimension d (finite) of the system increases, 
the asymptotic slope of Gy and Gp should approach the corresponding value for a Poisson 
point process in dimension d + 1. The results in Table [l] indicate that this claim closely 
holds for d = 3 and 4, meaning that the convergence of processes is relatively quick with 
respect to increasing dimension. Based on the analysis in p], we therefore expect this trend 
to continue for higher dimensions. 



d Gv 

1 7rV2 (exact) 

2 2.499 ±0.015 

3 1.680 ±0.025 

4 1.323 ±0.049 



TABLE I: Large-r slopes of Gy for each dimension. The d = 1 slope is taken from the asymptotic 



expansion in (57). Given errors are estimated based on the approximate error for d = 1. 



3. Voronoi statistics of the Fermi sphere point process for d = 2 

To demonstrate the utility of the HKPV Algorithm in statistically characterizing a point 
process, we have also included statistics for the Voronoi tessellation of the d = 2 Fermi sphere 
point process in Table [11} Specifically, we provide results for the probability distribution of 
the number of cell sides p„ and the average area of an n-sided cell (v4„). Similar results have 
been reported in the literature for Voronoi tessellations of Poisson point processes and 
determinantal point processes generated from the eigenvalues of complex random matrices 
[38] : we also provide the comparison in Table [ll[ Visual representations of the data are shown 
in Figure [TTj 

The topology of the plane enforces the constraints that (n) = 6 and (A) = 1/p (= 1 at 
unit density) for any point process, where n is the number of cell sides and A is the area of 
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TABLE II: Voronoi statistics for several point processes with d = 2. FPP = Fermi sphere point 
process; PPP = Poisson point process; CRM = complex random matrix. Results for the PPP and 
CRM are from [3^. The systems have been normalized to unit number density (p = 1). 




n 



FIG. 11: Left: Distribution p„, of the number of sides n of Voronoi cells for the Fermi sphere point 
process (FPP), Poisson point process (PPP), and eigenvalues of a complex random matrix (CRM). 
Right: Expectation value of the area of an n-sided Voronoi cell {An) for the FPP, PPP, and CRM. 

a cell. We notice that the distribution p„ is more sharply peaked for the Fermi sphere point 
process than in the Poisson point process, which is a consequence of the effective repulsion 
among the particles. With regard to the average areas of cells, is appears that Fermi-sphere 
cells with smaller n have larger areas than Poisson cells, again likely due to the repulsion of 
the points; however, Poisson cells with a greater number of sides tend to have larger areas 
than Fermi-sphere cells, a result which can be attributed to the more even distribution of 
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points in the Fermi-sphere process through space, which is related to the hyperuniformity of 



the point process. Figure 12 shows a typical Voronoi tessellation for the Fermi-sphere point 
process compared to the equivalent tessellation for a Poisson point process. We immediately 
notice that the determinantal point process tends to avoid clustering of particles, resulting in 
a narrower distribution of cell sizes within the tessellation; such clustering is not precluded 
in the Poisson tessellation, resulting in isolated regions of small (or large) cells. 



FIG. 12: Left: Voronoi tessellation of the d = 2 Fermi sphere point process at number density 
p = l. 

Right: Voronoi tessellation of a d = 2 Poisson point process at number density p = 1. 
Both: Tessellations are performed with periodic boundary conditions using N = 109 points. 

In order to rationalize these properties, we utilize the hyperuniformity (superhomogene- 
ity) of the Fermi sphere point process. Voronoi tessellations of hyperuniform point processes 
share several unique characteristics which distinguish them from general point processes. 
For example, Gabrielli and Torquato [IT] have provided the following summation rule, which 
holds for all hyperuniform point processes in any dimension: 

( E / = E = 0' (78) 

\ j=l / j=-oo 

where V is the system volume, N[S) is the number of points in a large subset S of V, 
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Wi = Vi — 1/p, Vi is the size of Voronoi cell i, and Cij = {wiWj) defines the correlation 
matrix between the sizes of different Voronoi cells. We note that this rule is essentially a 
discretization of the condition that S{0) = for a hyperuniform point process, meaning that 
infinite-wavelength number fluctuations vanish within the system; therefore, the result in 



(78) is unique to tessellations of hyperuniform point processes. Additionally, Gabrielli and 
Torquato have shown that arbitrarily large Voronoi cells or cavities are permitted in hyper- 
uniform point processes despite the fact that these processes possess the slowest growth of 
the local-density fluctuation with R (the size of the window) [T7j. We particularly emphasize 
the result that the probability distribution of the void regions must decay faster in R than 
the equivalent distribution for any non-hyperuniform process. 

We show elsewhere [2] that the structure factor in any dimension d for the Fermi sphere 
point process has the following nonanalytic behavior at the origin: 

S{k) r^k {k^ 0), (79) 

and the large- number variance is controlled by: 

a\R) ^ R'^-HniR). (80) 

The unusual asymptotic scaling a'^{R)/R'^~^ = ln(i?) for the Fermi sphere point process has 
also been observed in three-dimensional maximally random jammed sphere packings [39] , 
which can be viewed as prototypical glasses since they are both perfectly rigid mechanically 
and maximally disordered. 

The peaking phenomenon observed in the Voronoi statistics of the Fermi sphere point 
process therefore reflects the fact that the probability of observing large Voronoi cells must 
be less than the corresponding probability for the Poisson point process, which is not hyper- 
uniform. The more even distribution of the Voronoi cells through space in the Fermi sphere 
point process prevents the probability distribution of the cell sizes from decaying more slowly 
than the corresponding distribution for the Poisson point process, where clustering of the 
points increases the likelihood of observing both smaller and larger Voronoi cells. 

The comparison between the Voronoi statistics of the Fermi sphere point process and 



the Ginibre ensemble in Figure [TT] highlights the similarities between the two determinantal 
point processes. Namely, the distributions pn for each system are sharply peaked around 
n = 6 and narrower than the corresponding result for the Poisson point process. However, 
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notable differences between the statistics are also apparent. The distribution j9„ for the 
Fermi sphere point process is more sharply peaked than the corresponding result for the 
Ginibre ensemble. The larger probability in the Ginibre ensemble of observing cells with a 
fewer or larger number of sides n is directly related to the correlations among the particles 
in the system. 



C. Comparison of results across dimensions 



In order to compare statistical quantities across dimensions, it is generally preferable to 
enforce a fixed mean nearest-neighbor separation A since this quantity determines the length 
scale of the system. This constraint is easily obtained via a rescaling of the density according 
to the relation: 

Hp) = A(i) l^-j , (81) 



where A(l) denotes the mean nearest-neighbor separation at unit density. Equation (81) 
easily follows from the scaling of the density p with the size of the system. Of particular 
interest are the values of A(l) for each dimension and p(l), the number density at which 
the system has unit mean nearest-neighbor separation. These quantities may be read from 
Table [ml 



It is not difficult to show, using (81), that p(l) = A(l) . We note that for sufficiently 
large p, the mean nearest-neighbor separation increases with the dimension of the system; 
however, the opposite trend is observed for small p. For intermediate values of the density, 
the trend becomes less discernable. At unit density, we observe that A(l) decreases between 
d = 1 and d = 2 but then increases again for d >2; indeed, we measure this trend directly in 



Table III Estimates for A(l), which are developed elsewhere [2J, suggest that A(l) continues 
to increase with respect to increasing dimension; if this result is true, then we therefore 



expect that as d — > oo, p(l) oo. From the definitions of g2 and kp in (75), one can show 
that g2{r) = g^\X{l)r], where g^^ is the form of the pair correlation function at unit density. 
Therefore, as A(l) increases, the curve representing g2 shifts to the left, implying that for 
large dimensions g2 is approximately given by unity for all r, and the system is uncorrelated. 
This behavior is a direct consequence of enforcing a fixed mean nearest-neighbor separation 
on the system as opposed to a fixed density. 



After appropriate rescaling, we compare the results for Gp and Gy in Figure 14 The 
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1 0.725728 0.725728 

2 0.649823 0.422270 

3 0.654511 0.280382 

4 0.679561 0.213262 



TABLE III: Values of A(l) and p(l) for each dimension. 

results strongly suggest that Gy{r) — 1 as (i — > oo, which is in agreement with the con- 
clusions drawn from the analysis above. We also notice that both Gp and Gy decrease in 
slope as the dimension of the system increases; thus, if Gp and Gy possess the same r ^ oo 
asymptotic slope, then it must be true that Gp saturates at unity for large r in the limit 
d —>■ oo. This behavior is surprising in the context of our description of g2 above. The fact 
that (72 1 for large d indicates a decorrelation of the system for higher-dimensions, leading 
us to expect Poisson-like behavior in the system as conjectured in [2]. The behavior of Gy 
corroborates this notion as does the convergence of Gp and Gy for large d. However, our un- 
derstanding of Hp and Ep from the discussion above along with the bounds from [2], which 
sharpen with increasing dimension at fixed A, suggest instead that Hp Hp = 6{r — 1) 
and Ep Ep = 0(1 — r) for large d, where S{x) is the Dirac delta function, G(x) is the 
Heaviside step function, and Hp and Ep are effective generalized functions. As shown in 
[2], the only functional form for Gp that agrees with these conclusions and the observed 



behavior in Figure 14 is Gp G*p = 0(r — 1) as — oo for fixed mean nearest-neighbor 
separation. 

We rationalize these observations by noting that the effective hard core of the fermionic 
system as described in [2] has been encoded in the functional form of Gp due to the constraint 
of fixed mean nearest-neighbor separation. It is this constraint which produces the limiting 
forms of Hp and Ep for high dimensions, meaning that the environment around any given 
particle greatly resembles a saturated system of hard spheres. However, the scaling of g2 
with A(l) mentioned above means also that the particles only see large-r correlations from 
the corresponding form of g2 at unit density, resulting in Poisson-like behavior for this 
function [2j, which is then translated into the value of unity for Gy in high dimensions. In 



36 




FIG. 13: Left: Hp(s) for the Fermi sphere point process at unit mean nearest-neighbor separation 
A for d = 1,2, 3, 4. Right: Ep{s) for the Fermi sphere point process at unit mean nearest-neighbor 
separation A for d = 1,2, 3, 4. 





FIG. 14: Comparison of Gp {left) and Gy (right) across dimensions at unit mean nearest-neighbor 



separation A. Results are from numerical calculations using (19) 



other words, the particle quantities contain the information about the effective hard core 
under the constraint of fixed mean nearest-neighbor separation, but the void quantities are 
Poisson-hke to account both for the scahng of g2 and the small- and large-r contraints shown 
numerically in Section IV. B that must be enforced regardless of how the infinite-dimensional 
limit is taken. 



Figure 15 shows the distributions of the extremum nearest-neighbor distances at fixed 
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FIG. 15: Distributions of the maximum (M*(s); left) and minimum (M=k(s); right) nearest-neighbor 
distances across dimensions at unit mean nearest-neighbor separation A. Results are simulated 
using the HKPV Algorithm. 



mean nearest-neighbor separation based on calculations from configurations generated with 
the HKPV Algorithm. We have been unable to write these quantities in determinantal 
form amiable to numerical calculation, and therefore the HKPV Algorithm is an attractive 
means through which to study these quantities. We note that the maximum and minimum 
nearest-neighbor spacings appear to converge to a value of unity as the dimension of the 
system increases; this behavior is expected in the context of the discussion for Hp above. 



The convergence of these quantities is more easily seen in Figure 16 we have also included 
the values of A(l) for reference, but there is strong evidence to suggest that the limiting 
value of the extremum quantities for large d is unity. 

V. DETERMINANTAL PROCESSES IN CURVED SPACES 

In this last section, we present an example of how the HKPB algorithm is not limited to 
point processes in Euclidean spaces described above. With an appropriate choice of the basis 
functions 0n it can in principle be adapted to simulate point processes on other domains 
and topologies. Of particular interest in this regard is the generation of point processes 



on a curved space, like the two-sphere 5* in Figure 17 Here, we consider the spherical 



harmonics as basis functions for a spherical geometry; 0„ = Yi^rn{d,(p) are a basis for the 
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— • Minimum nearest-neighbor distance 
— • Maximum nearest-neiglibor distance 
• Unit mean nearest-neighbor distance 
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FIG. 16: Average maxima and minima nearest-neighbor distances with standard deviations across 
dimensions at unit mean nearest-neighbor separation; results are obtained using the HKPV Algo- 
rithm. Also included for reference is the unit mean nearest-neighbor separation, which is fixed for 
each dimension. 




FIG. 17: Configuration of 37 points on the unit sphere using the HKPB algorithm with the spherical 
harmonics as basis functions. 




39 



square-integrable functions on the two-sphere 5*^. Since m = —I, so for any / there are 
2/ + 1 different values of m, we decided to choose the lowest (A^ — l)/2 values of I and all 
the corresponding m's. Once these functions have been chosen, the algorithm provides a 
relatively simple means to generate the point process. We have not embarked in an extensive 
analysis of the statistical properties of this process as we leave that for future work. We note, 
however, from previous observations that a short-distance effective interaction among the 
points is logarithmic and repulsive, and we expect a fluid-like configuration on the surface 
of the sphere. Also, for — > oo at fixed sphere radius it is not difficult to conjecture that 
the nearest neighbor functions will tend to those we already discussed for the Fermi-sphere 
process on torus. On the other hand, for finite A^ this problem could be relevant to problem 
of packing of spheres in non- Euclidean geometries. This is a promising direction for future 
research. 



VI. CONCLUDING REMARKS 



Our focus in this paper has been on characterizing the statistical properties of high- 
dimensional determinantal point processes through both numerical calculations and algo- 
rithmic generations of point configurations. We first compared the results for n-particles 
distribution functions and nearest-neighbor function obtained by the two methods to cross- 
check consistency and accuracy. We then proceeded using both methods to elucidate the 
small- and large-r behaviors of the nearest-neighbor distribution functions and the extrema 
statistics in dimensions one to four. Our results strongly suggest that both Gp and Gy are 
linear for sufficiently large r, and we obtain numerical estimates of the common slope in this 
limit. This behavior is to be contrasted with the equivalent forms of Gp and Gy for equi- 
librium systems of hard spheres and for Poisson point processes. It is known for the former 
system that both functions saturate for sufficiently large r while Gp{r) = Gv{r) = 1 for all 
r in the latter process [4J. The linearity of Gp and Gy in the determinantal case is thus 
unique in the context of general point processes. We have also shown, in accordance with 
[2], that in the limit as ci oo, both Gp and Gy must saturate at unity in accordance with 
the behavior of the aforementioned bounds; again, this claim is supported by the numerical 
evidence we have presented here. Also, as the dimension d grows, we observed that the 
functions Hp and Hy become concentrated around their maximum as do the distributions 
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of extrema of nearest-neighbor distances and M* . 

By using the HKPV algorithm to generate configurations of points we have shown that 
the determinantal nature of the n-particle probabihty density has a significant effect on 
the Voronoi cell statistics of the Fermi-sphere point process in two dimensions. Namely, 
the probability distribution of cell-sides is more peaked around n = Q (hexagons) than the 
corresponding distribution for either Poisson point process or the Ginibre ensemble [38] (the 
distribution of complex eigenvalues of random complex matrices). The effective separation 
of the points, resulting in a sharper peak in the distribution of the number of sides of the 
Voronoi cells, is closely related to the hyperuniformity of the system. 

Finally, to show how the algorithm can be used for generating determinantal processes 
on curved spaces, we have presented an example of determinantal point process on the 
two- sphere. 
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